## Replication File for:
## ``Concentrated Burdens: How Self-Interest and Partisanship Shape Opinion on Opioid Treatment Policy"
## Justin de Benedictis-Kessner and Michael Hankinson


## ---------- ##
#### Notes: ####
## ---------- ##

# This is the first script file needed to reproduce the results from the paper
# Run this script first to produce recoded survey dataset that subsequent scripts need
# Other scripts will produce the tables and figures that represent the main results


## --------------------- ##
#### Preliminary stuff ####
## --------------------- ##

rm(list=ls())
library(foreign)
library(car)
library(descr)
library(xtable)
library(ggplot2)
library(reshape2)
library(stargazer)
library(tidyverse)


## ----------------- ##
#### Read in data: ####
## ----------------- ##
a1 <- read_csv("s578_client_20180906.csv") # Data collected August 2018 via NORC AmeriSpeak Panel

# add county level racial and segregation data
county <- read_csv("zips_county_race.csv")
county2 <- read_csv("R11814798_SL050.csv")
county$pop_density <- county2$SE_T002_002[match(county$county_id,county2$Geo_FIPS)]

a <- merge(a1, county, by.x="ZIP_CODE", by.y="zip5", all.x=T)


# ------------------ #
#### Data Recodes ####
# ------------------ #

## add sample death rate binary variable
# mean(a$death_rate, na.rm=T)
a$death_rate_bi <- ifelse(a$death_rate>8,1,0)

# definitions
a$condition <- ifelse(a$P_BASED=="2","ability", "need")
a$above_death_rate <- ifelse(a$P_ZIPREF=="1",1,0)
a$above_income <- ifelse(a$P_INCREF=="1",1,0)

a$ability_support <- ifelse(a$condition=="ability", a$Q1, NA)
a$needs_support <- ifelse(a$condition=="need", a$Q1, NA)

# convert to numeric unit scale. '98' as skipped
a$ability_support <- car::recode(as.numeric(a$ability_support),"1=1; 2=.75; 3=.5; 4=.25; 5=0; 98=NA")
a$needs_support <- car::recode(as.numeric(a$needs_support),"1=1; 2=.75; 3=.5; 4=.25; 5=0; 98=NA")

### support as binary variable
# 0 = somewhat/strongly oppose, 1 = strongly/somewhat support + indifferent
a$ability_support_bi <- as.numeric(a$ability_support>0.5)
a$needs_support_bi <- as.numeric(a$needs_support>0.5)

# gender:
a$female <- as.numeric(a$GENDER==2)

# PID. Some negative numbers. -1 Dem, 0 Independent, 1 Republican
# hist(a$P_PARTYID7)
a$P_PARTYID7 <- ifelse(a$P_PARTYID7>0, a$P_PARTYID7, NA) #assume -1 means did not answer
a$pid <- ifelse(a$P_PARTYID7 < 4,-1, ifelse(a$P_PARTYID7 > 4, 1, 0))
a$dem <- as.numeric(a$pid==-1)
a$rep <- as.numeric(a$pid==1)

# Ideology, Some negative number and 8 = Haven't thought about it. -1 Liberal, 0 Moderate, 1 Conservative
# hist(a$P_IDEO)
a$P_IDEO <- ifelse(a$P_IDEO<0,NA,ifelse(a$P_IDEO>7, NA, a$P_IDEO) )
a$ideo <- ifelse(a$P_IDEO<4,-1, ifelse(a$P_IDEO>4, 1, 0))
a$conservative <- as.numeric(a$ideo==1)
a$liberal <- as.numeric(a$ideo==-1)

# income:
# hist(a$INCOME)
a$income_gt50k <- as.numeric(a$INCOME>9)

# education:
# hist(a$EDUC4)
a$college <- as.numeric(a$EDUC>=11) # college degree of some kind

# marital status:
# table(a$MARITAL)
a$married <- as.numeric(a$MARITAL==1)

# race:
a$white <- as.numeric(a$RACETHNICITY==1) 
a$black <- as.numeric(a$RACETHNICITY==2) 
a$hispanic <- as.numeric(a$RACETHNICITY==4)

# home ownership
# table(a$HOUSING)
a$own <- ifelse(a$HOUSING=="1", 1, ifelse(a$HOUSING=="2", 0, NA))

# population density - binned
a$pop_density_bin <- cut(a$pop_density,quantile(a$pop_density,c(0,0.25,0.5,0.75,1),na.rm=T))
a$pop_density_bin_num <- as.numeric(a$pop_density_bin)
# hist(a$pop_density_bin_num)

# create combined support variable preference from both conditions:
a$support <- ifelse(a$condition=="ability", a$ability_support, a$needs_support)
a$support_bi <- ifelse(a$condition=="ability", a$ability_support_bi, a$needs_support_bi)
# mean(a$support_bi, na.rm=T)

## Personal exposure
# table(a$Q3)
nrow(subset(a, Q3==1))/nrow(a) # .03 personally addicted
1-(nrow(subset(a, Q3==5))/nrow(a)) # .54 know someone who has been addicted

# convert to numeric unit scale, tossing out 98s (refusal to answer?)
a$personal <- car::recode(as.numeric(a$Q3),"1=1; 2=.75; 3=.5; 4=.25; 5=0; 98=NA")
# table(a$personal)
# exposure as binary variable
# 0 = no knowledge of anyone, 1 = me, family, friend, acquaintance
a$personal_bi <- as.numeric(a$Q3<5)

### NIMBYism treatment condition:
a$distance <- ifelse(a$P_DISTANCE=="2","far", "near")

a$near_support<-ifelse(a$P_DISTANCE=="1", a$Q2, NA)
a$far_support<-ifelse(a$P_DISTANCE=="2", a$Q2, NA)

# convert to numeric unit scale:
a$near_support <- car::recode(as.numeric(a$near_support),"1=1; 2=.75; 3=.5; 4=.25; 5=0; 98=NA")
a$far_support <- car::recode(as.numeric(a$far_support),"1=1; 2=.75; 3=.5; 4=.25; 5=0; 98=NA")

### support as binary variable
# 0 = somewhat/strongly oppose, 1 = strongly/somewhat support + indifferent
a$near_support_bi <- as.numeric(a$near_support>0.5)
a$far_support_bi <- as.numeric(a$far_support>0.5)

# create combined support variable preference from both conditions:
a$dist_support <- ifelse(a$distance=="near", a$near_support, a$far_support)
a$dist_support_bi <- ifelse(a$distance=="near", a$near_support_bi, a$far_support_bi)

a$state_fips <- as.numeric(gsub("(.{1,2})\\d{3}","\\1",a$county_id))

# region variables:
northeast <- c(09,23,25,33,44,50,34,36,42)
south <- c(10,11,12,13,24,37,45,51,54,01,21,28,47,05,22,40,48)
west <- c(04,08,16,30,32,35,49,56,02,06,15,41,53)
midwest <- c(18,17,26,39,55,19,31,20,38,27,46,29)

a$northeast <- as.numeric(a$state_fips %in% northeast)
a$south <- as.numeric(a$state_fips %in% south)
a$west <- as.numeric(a$state_fips %in% west)
a$midwest <- as.numeric(a$state_fips %in% midwest)

a$region <- recode(a$REGION4,`1` = "Northeast",
									 `2` = "Midwest",
									 `3` = "South",
									 `4` = "West"
)

a <- a %>%
	select(-ZIP_CODE) # removing for confidentiality per NORC DUA

write_csv(a,"norc_data_recoded.csv")
